### Let's make some pretty pictures
### Need to produce a 4-paneled figure for different topographic features (D2C, DEM, BVG, Aspect, Slope, etc...)
### Panels will be freq distributions (histograms) for the whole region, within rainforest, sampled sites, and a comparison between sampled sites and rainforest

#list the libraries needed
necessary=c("adehabitat","SDMTools","rgdal","sp","Hmisc")

#check if library is installed
installed = necessary %in% installed.packages()

#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)

#load the libraries
for (lib in necessary) library(lib,character.only=T)

# Establish directories

out.dir = '/home1/99/jc152199/TopoGapAnalysis/images/'
in.dir = '/home1/99/jc152199/TopoGapAnalysis/ASCIIFiles/'
setwd(in.dir)

# Read in ASCII files

slope.asc=read.asc(paste(in.dir,"slope80WTplusbuffer.asc",sep=""))
aspect.asc=read.asc(paste(in.dir,"aspect80WTplusbuffer.asc",sep=""))
coastdist.asc=read.asc(paste(in.dir,"coastdist80WTplusbuffer.asc",sep=""))
dem.asc=read.asc(paste(in.dir,"dem80WTplusbuffer.asc",sep=""))
BVG.asc=read.asc(paste(in.dir,"BVG2M80WTplusbuffer.asc",sep=""))

# Create a matrix of unique row/column positions for the whole region

base.pos = as.data.frame(which(is.finite(BVG.asc), arr.ind = T))
base.pos$east = getXYcoords(BVG.asc)$x[base.pos$row]
base.pos$north = getXYcoords(BVG.asc)$y[base.pos$col]

# Bind row/column info and ASCII information into a single data frame
# base.pos will now be the object upon which regional frequency distributions of topography are based

base.pos$BVG = BVG.asc[cbind(base.pos$row,base.pos$col)]
base.pos$slope = slope.asc[cbind(base.pos$row,base.pos$col)]
base.pos$aspect = aspect.asc[cbind(base.pos$row,base.pos$col)]
base.pos$dem = dem.asc[cbind(base.pos$row,base.pos$col)]
base.pos$coastdist = coastdist.asc[cbind(base.pos$row,base.pos$col)]

# Read in survey loc info

surveylocs = read.csv(paste('/home1/99/jc152199/TopoGapAnalysis/cjs_Rsurveylocs.csv',sep=""),header=T)
survey.pnts = cbind(surveylocs$east,surveylocs$north)

# Use extract.data command to get info from the ASCII files for all survey locations

survey.slope = extract.data(survey.pnts, slope.asc)
survey.coastdist = extract.data(survey.pnts, coastdist.asc)
survey.aspect = extract.data(survey.pnts, aspect.asc)
survey.dem = extract.data(survey.pnts, dem.asc)
survey.BVG = extract.data(survey.pnts, BVG.asc)
survey.BVGRFonly = extract.data(survey.pnts, BVGRFonly.asc)

survey.topo = data.frame(east=survey.pnts[,1], north=survey.pnts[,2], coastdist = survey.coastdist, slope = survey.slope, aspect = survey.aspect, dem = survey.dem)


# Identify all positions within the region that are within the veg types of interest (1-8, all rainforest types and Wet Sclerophyll

veg.pos = base.pos[which(base.pos$BVG<=8),]

veg.pnts = cbind(veg.pos$east, veg.pos$north)

# Use extract.data command to get info from the ASCII files for all rainforest locations
# This data will be used for the frequency distribution of topographic features within rainforest

veg.slope = extract.data(veg.pnts, slope.asc)
veg.coastdist = extract.data(veg.pnts, coastdist.asc)
veg.aspect = extract.data(veg.pnts, aspect.asc)
veg.dem = extract.data(veg.pnts, dem.asc)
veg.BVG = extract.data(veg.pnts, BVG.asc)

# Produce histograms of slope for whole region

png(filename=paste(out.dir,"slope4panel.png",sep=""), units="cm", res=100, height=40, width=120, bg="white")

par(mfrow=c(1,3), cex=2)

hist(slope.asc, freq=F, ylim=c(0,.2), xlim=c(0,round(max(slope.asc,na.rm=T),0)+3), breaks=seq(round(min(slope.asc, na.rm=T),0),round(max(slope.asc, na.rm=T),0)+3,by=3))
hist(veg.slope, freq=F, ylim=c(0,.2), xlim=c(0,round(max(slope.asc,na.rm=T),0)+3), breaks=seq(round(min(slope.asc, na.rm=T),0),round(max(slope.asc, na.rm=T),0)+3,by=3))
hist(survey.slope, freq=F, ylim=c(0,.2), xlim=c(0,round(max(slope.asc,na.rm=T),0)+3), breaks=seq(round(min(slope.asc, na.rm=T),0),round(max(slope.asc, na.rm=T),0)+3,by=3))

dev.off()

# Produce histograms of aspect for whole region

png(filename=paste(out.dir,"aspect4panel.png",sep=""), units="cm", res=100, height=40, width=120, bg="white")

par(mfrow=c(1,3), cex=2)

hist(aspect.asc, freq=F, ylim=c(0,.01), xlim=c(0,360), breaks=seq(round(min(aspect.asc, na.rm=T),0),round(max(aspect.asc, na.rm=T),0)+15,by=15))
hist(veg.aspect, freq=F, ylim=c(0,.01), xlim=c(0,360), breaks=seq(round(min(aspect.asc, na.rm=T),0),round(max(aspect.asc, na.rm=T),0)+15,by=15))
hist(survey.aspect, freq=F, ylim=c(0,.01), xlim=c(0,360), breaks=seq(round(min(aspect.asc, na.rm=T),0),round(max(aspect.asc, na.rm=T),0)+15,by=15))

dev.off()

# Produce histograms of coastdist for whole region

png(filename=paste(out.dir,"coastdist4panel.png",sep=""), units="cm", res=100, height=40, width=120, bg="white")

par(mfrow=c(1,3), cex=2)

hist(coastdist.asc, freq=F, ylim=c(0,.00004), xlim=c(0,90000), breaks=seq(round(min(coastdist.asc, na.rm=T),0),round(max(coastdist.asc, na.rm=T),0)+5000,by=5000))
hist(veg.coastdist, freq=F, ylim=c(0,.00004), xlim=c(0,90000), breaks=seq(round(min(coastdist.asc, na.rm=T),0),round(max(coastdist.asc, na.rm=T),0)+5000,by=5000))
hist(survey.coastdist, freq=F, ylim=c(0,.00004), xlim=c(0,90000), breaks=seq(round(min(coastdist.asc, na.rm=T),0),round(max(coastdist.asc, na.rm=T),0)+5000,by=5000))

dev.off()

# Produce histograms of dem for whole region

png(filename=paste(out.dir,"dem4panel.png",sep=""), units="cm", res=100, height=40, width=120, bg="white")

par(mfrow=c(1,3), cex=2)

hist(dem.asc, freq=F, ylim=c(0,.003), xlim=c(0,1600), breaks=seq(round(min(dem.asc, na.rm=T),0),round(max(dem.asc, na.rm=T),0),by=100))
hist(veg.dem, freq=F, ylim=c(0,.003), xlim=c(0,1600), breaks=seq(round(min(dem.asc, na.rm=T),0),round(max(dem.asc, na.rm=T),0),by=100))
hist(survey.dem, freq=F, ylim=c(0,.003), xlim=c(0,1600), breaks=seq(round(min(dem.asc, na.rm=T),0),round(max(dem.asc, na.rm=T),0),by=100))

dev.off()


# Mess around with aggregate command to create a chart of the differences between topographic frequency distributions within rainforest and at sample sites

a= data.frame(east=survey.pnts[,1], north=survey.pnts[,2], slope=round(survey.slope,0))
b=aggregate(a$east, by=list(a$slope), FUN=length)
d=sum(b$x,na.rm=T)
b$perc = NULL
b$perc = b$x/d

e = data.frame(east=veg.pnts[,1], north=veg.pnts[,2], slope = round(veg.slope,0))
f = aggregate(e$east, by=list(e$slope), FUN=length)
g = sum(f$x, na.rm=T)
f$perc = NULL
f$perc = f$x/g

# Turn the BVG ASCII into a binary matrix (1 represents RF 0 represents all other veg types)

BVGbinary.asc = BVG.asc
BVGbinary.asc[which(BVG.asc>=9)]=0  # Other veg types
BVGbinary.asc[which(BVG.asc<=8)]=1  # Rainforest and Wet Sclerophyll

write.asc(x=BVGbinary.asc, file=paste(out.dir,'BVGbinary.asc',sep=""))

# Image of resulting ASCII

png(filename=paste(out.dir,"BVGbinary.png",sep=""), units="cm", res=100, height=80, width=80, bg="white")
image(BVGbinary.asc)
dev.off()

# Run ConnCompLabel to identify individual RF patches within the landscape

BVGpatches = ConnCompLabel(BVGbinary.asc)

# BVGpatches[which(BVGpatches==0)]=NA

write.asc(x=BVGpatches, file=paste(out.dir,"patches.asc",sep=""))

# Link each row/column from base.pos with the unique id of it's patch from BVGpatches

base.pos$patchid = BVGpatches[cbind(base.pos$row,base.pos$col)]

# Run PatchStat will holds information for each patch in the landscape based on it's unique ID
# Problem ### PatchStat computes '0' from BVGBinary (non-rainforest) as a patch, and consequently fucks up and computes the largest RF Patch as the non-rainforest

z = PatchStat(BVGpatches)
z[1,]=NA # This command takes the first row of patch stats (the non-rainforest areas, and makes all their patch stats NA)

# Merge patch info from 'z' with base.pos based on unique patch id

y = merge(base.pos,z,by.x='patchid',by.y='patchID',all.x=F)
y$logarea = log(y[,10])

# Create a new blank ASCII

RFPatchSize.asc = BVG.asc                                                                                              

RFPatchSize.asc[,] = NA

# Populate blank ASCII with patch data from 'y'

RFPatchSize.asc[cbind(y$row,y$col)] = y[,10]

write.asc(x=RFPatchSize.asc, file=paste(out.dir,'RFPatchSize.asc',sep=""))

png(filename=paste(out.dir,"RFpatchsize.png",sep=""), units="cm", res=100, height=40, width=40, bg="white")

#par(mfrow=c(1,3), cex=2)

hist(y$logarea)

dev.off()

# Get survey data of RF Patch Sizes

a = extract.data(survey.pnts, RFPatchSize.asc)
b = log(a)

png(filename=paste(out.dir,"RFpatchsizeSurveysOnly.png",sep=""), units="cm", res=100, height=40, width=40, bg="white")

#par(mfrow=c(1,3), cex=2)

hist(b)

dev.off()

# Fuck around with MDS & PCA

pca = prcomp(survey.topo)
pca2 = prcomp(survey.topo[,3:6])

distmat = dist(survey.topo, method="euclidean", diag=T, upper=T)
zz = cmdscale(distmat,k=2,eig=T, x.ret=T)

png(filename=paste(out.dir,"mdstest.png",sep=""), units="cm", res=100, height=40, width=40, bg="white")

plot(zz$points[,1], zz$points[,2], xlab="Dimension 1", ylab="Dimension 2",main ="MDS")    #put up a graphics window

dev.off()

png(filename=paste(out.dir,"pcabiplot2.png",sep=""), units="cm", res=100, height=40, width=40, bg="white")

biplot(prcomp(survey.topo[,3:6]))    #put up a graphics window

dev.off()

vegdistmat = vegdist(survey.topo, method='euclidean')



